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Abstract 

A stochastic representation for the solutions of the Poisson-Vlasov 
equation, with several charged species, is obtained. The representation 
involves both an exponential and a branching process and it provides 
an intuitive characterization of the nature of the solutions and its 
fluctuations. Here, the stochastic representation is also proposed as a 
tool for the numerical evaluation of the solutions 

1 Introduction 

It is well known that the solutions of linear elliptic and parabolic equations, 
both with Cauchy and Dirichlet boundary conditions, have a probabilistic 
interpretation. This is a very classical field which may be traced back to 
the work of Courant, Friedrichs and Lewy [1] in the 20's. In spite of the 
pioneering work of McKean [2|, the question of whether useful probabilistic 
representations could also be found for a large class of nonlinear equations 
remained an essentially open problem for many years. It was only in the 90's 
that, with the work of Dynkin^ such a theory started to take shape. 
For nonlinear diffusion processes, the branching exit Markov systems, that 
is, processes that involve both diffusion and branching, seem to play the same 
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role as Brownian motion in the linear equations. However the theory is still 
limited to some classes of nonlinearities and there is much room for further 
mathematical improvement. 

Another field, where considerable recent advances were achieved, was the 
probabilistic representation of the Fourier transformed Navier-Stokes equa- 
tion, first with the work of LeJan and Sznitman[S], later followed by extensive 
developments of the Oregon schoolfB] [Z| [H]. In all cases the stochastic rep- 
resentation defines a process fort which the mean values of some functionals 
coincide with the solution of the deterministic equation. 

Stochastic representations, in addition to its intrinsic mathematical rele- 
vance, have several practical implications: 

(i) They provide an intuitive characterization of the equation solutions; 

(ii) By the study of exit times from a domain they sometimes provide 
access to quantities that cannot be obtained by perturbative methods [9j 

(iii) They provide a calculation tool which may replace, for example, the 
need for very fine integration grids at high Reynolds numbers; 

(iv) By associating a stochastic process to the solutions of the equation, 
they may also provide an intrinsic characterization of the nature of the fluctu- 
ations associated to the physical system. In some cases the stochastic process 
is essentially unique, in others there is a class of processes with means leading 
to the same solution. 

In [TU] a stochastic representation has been obtained for the solutions of 
the Fourier-transformed Poisson-Vlasov equation in 3 dimensions for parti- 
cles of one charge species on an arbitrary background. Here this result is 
generalized for the case of several charged species. As before the represen- 
tation involves both an exponential and a branching process, the solution 
being obtained from the expectation value of a multiplicative functional over 
backwards in time realizations of the process. 

The backwards in time realization of the process turns out to be appropri- 
ate for (parallelizable) numerical evaluation of the solutions and the Fourier 
representation adequate to obtain information on the small scale behaviour. 

2 The stochastic representations 

Consider a multi-species Poisson-Vlasov equation in 3+1 space-time dimen- 
sions 
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{i = 1,2), with 



Passing to the Fourier transform 
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with rj = i^x^v^ and ^ = ^^1)^2^ = (^1)^2)5 one obtains 
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Changing variables to 

where 7(|^2|) is a positive continuous function satisfying 
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with ^1= il^. Eq.Q written in integral form, is 
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For convenience, a stochastic representation is going to be written for the 
following function 
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with A a constant and h{^i) a positive function to be specified later on. The 
integral equation for x (^ii "T") is 
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Eq.Q has a stochastic interpretation as an exponential process (with a 
time shift in the second variable) plus a branching process. P {^.i, Ci) d^Ci 
is the probability that, given a mode, one obtains a (^i — ^1,^1) branch- 
ing with ^[ in the volume (^1,^1 + d^^'i)- xi£,ii^2:'^) is computed from the 
expectation value of a multiplicative functional associated to the processes. 
Convergence of the multiplicative functional hinges on the fulfilling of the 
following conditions : 



(A) ^ < 1 

(B) (|«;r'fc*/j)«i) <A({i) 

Condition (B) is satisfied, for example, for 



and 
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Indeed computing ^ h* h one obtains 
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Then ^ h* {^i) is bounded by a constant for all |^i|, and choos- 

ing c sufficiently small, condition (B) is satisfied. 

Once h {^i) consistent with (B) is found, condition (A) only puts restric- 
tions on the initial conditions. Now one constructs the stochastic process 

Because e"^"^ is the survival probability during time r of an exponential 
process with parameter A and Ae~^^(is the decay probability in the interval 
{s,s + ds), ^2) ■T") in Eq-® is obtained as the expectation value of a 

multiplicative functional for the following backward-in-time process, which 
we denote as process I : 

Starting at (^i,^2,t), a particle of species i lives for an exponentially 
distributed time s up to time t — s. At its death a coin Ig (probabilities |, ^) 
is tossed. If = two new particles of the same species as the original one 

are born at time r — s with Fourier modes — ^[,^2 + '^ ^(fgal) ) (^i' ^) 
with probability density If = 1 the two new particles are of 

different species. Each one of the newborn particles continues its backward- 
in-time evolution, following the same death and birth laws. When one of 
the particles of this tree reaches time zero it samples the initial condition. 
The multiplicative functional of the process is the product of the following 
contributions: 

- At each branching point where two particles are born , the coupling 
constant is 

V / niiX h{^i) 7(161) 

- When one particle reaches time zero and samples the initial condition 
the coupling is 

goi (6, 6) = — — (15) 
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The multiplicative functional is the product of all these couplings for each 
realization of the process X {^1,^2, ''"), this process being obtained as the limit 
of the following iterative process 
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Then, each Xi (Cij 61 ''') is the expectation value of the functional. 
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For example, for the realization in Fig.l the contribution to the multi- 
plicative functional is 



9ij sij gji [Ci- - S2j gu [ii^ii , r - 

xgoi (Cl - Ci', ^3) goi (^1 , 0) goj (^1,0) goi (^Ci - ^1 - ^1 , (18) 



and 



ki = /c + (r - Ti) 



^2 = A;i + (t2 - Ti) 

^3 = (7-3-Ti)Cl 

With the conditions (A) and (B), choosing 
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the absolute value of all coupling constants is bounded by one. The branch- 
ing process, being identical to a Galton- Watson process, terminates with 
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Figure 1: A sample path of the stochastic process I 



probabihty one and the number of inputs to the functional is finite (with 
probability one). With the bounds on the coupling constants, the multi- 
plicative functional is bounded by one in absolute value almost surely. 

Once a stochastic representation is obtained for x {Cii ^2, t), one also has, 
by ([8]), a stochastic representation for the solution of the Fourier-transformed 
Poisson-Vlasov equation and one obtains: 

Proposition 1. The process I, above described, provides a stochastic 
representation for the Fourier-transformed solutions of the Poisson-Vlasov 
equation Fi ((^1, ^2, t) for any arbitrary finite value of the arguments, provided 
the initial conditions at time zero satisfy the boundedness conditions (A). 

So far we have constructed a general process that provides a stochastic 
representation for the interacting Vlasov equation, not only for the Pois- 
son case, but also for more general situations with quadratic nonlinearities. 
However, because of the integrated nature of the Coulomb interaction, the 
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Poisson case is special in that there is also a representation by a simpler pro- 

cess. Looking at equation ([9]) one sees that because of the factor ,^2 ■ Ci ordy 
the trees where the mode Xj (^i)O)''" ^ ^) survives until time zero will con- 
tribute to the functional. That is, the only trees with non-zero contributions 
to the functional (fT7|) are the one-sided trees represented in Fig. 2. Therefore 
for the calculation of the solution one may replace Xj ('^i)0,r — s) by the 
initial condition computed at (^^[, (r — s) The process then becomes the 
following linear process with random couplings 

e 



the coupling constant at the branchings being 

(23) 

The functional representing the solution is the product of all branching 
coupling constants times one additional factor corresponding to the last non- 
branching mode. The result is the following 

Proposition 2. The linear process II, defined by / f^) and ( f^) also 

provides a stochastic representation of the solutions of the Poisson- Vlasov 
equation, the conditions on the kernels and initial conditions being given by 
(A), IHj and H^j. 



3 Stochastic representation and numerical codes 

The backwards-in-time probabilistic representations, obtained in Sect. 2, seem 
appropriate for the numerical evaluation of Fourier-localized solutions. Good 
statistics requires the average of the multiplicative functional over many re- 
alization trees. In the backwards in time realization one fixes a particular 
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Figure 2: A one-sided tree corresponding to process II 

mode at time r and generates as many trees as needed for that particular 
mode. Notice that by studying high Fourier modes one may obtain informa- 
tion about the small scale behaviour of the solution without having the need 
for a fine grid as it would be necessary in a real space numerical code. Each 
realization tree being independent of all the others, the probabilistic code is 
also appropriate for parallelization. 

We will not report, in this paper, extensive calculations using these rep- 
resentations and the corresponding codes. Nevertheless we list all the needed 
probability distributions needed to implement the method. 

For the construction of the sample trees and the calculation of the func- 
tional, the following probability densities are needed: 

- The probability of a ^i— mode branching into C,i and — modes 

with r ((^i) given by Eq. (IT^ . One notices that, for each this probability 
is only function of two variables, the \^[\ and the angle between and 
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Therefore defining 



^ = ^ (25) 

i+|e;r 



and changing the integration measure one obtains a 
- Probabihty density (z, cos ff) 

TT 1 

V {z, cos = —— (26) 

r(|6l) (^|^^|2 + i_2|^,|eos^yrrT) 

with z in the interval (0, 1) and cos^ in the interval (—1, 1). 

Because the inverse of the cumulative distribution functions have not a 
nice analytic form, we may use the reject method in the plane (2;,cos6') to 
simulate this probability distribution. For this, one needs 

P(^>cos^)^^ = ^^ (27) 



which is obtained for cos^ 6 = 1 and - = |^i|^ + 1. 

However because p {z , cos 9) ^^^^ is very large for large |^i| in a narrow 
region, it is more efficient to use 

- The integrated p (z) density, 

. . 2n 1 

with p (z)^^ = p (min (|6|~^ , l)), to choose z and then, once z is chosen, 
to use 

- The conditional probability density p(cos^|2;) 

p (cos 6*12;) = ^ — o (29) 



2(|6l' + i-2|ei|cos^yi^' 



with 



maxp (cos t/pj = ^ — ^ (cJUj 



2(l6r + i-2|6|^[^)" 



to choose cos 9. 
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Once z and cos^ are chosen, one computes = y ^ — 1 and chooses 

ip = 2'kRAND 

RAND being a random variable uniformly distributed in the interval (0,1) 
Then one obtains 



[sm U cos v?, sm U sm 99, cos ( 



and ^1 - 

Notice that only the amplitudes of the Fourier modes and |^i — 
are needed as inputs to compute the probabilities in the next branchings but 
the full vector is needed to compute ^2- Finally the lifetime r of each mode 
is obtained from 

A 

For the trees a standard indexation is used, each tree being a row vector 
of integer numbers with the number k at the position n meaning that that 
mode was born at the branching of the mode in the position k. 

The conditions (A), (fT2|) and (!20|) guarantee that all factors entering the 
multiplicative functional (ITTI) are bounded by one, implying that the func- 
tional itself is also bounded. This, together with the Galton- Watson nature 
of the branching, insures convergence of the expectation value. However, in 
practice, this leads to very small values of the functional and for large times 
(large trees) one may be faced with round-off inaccuracies in the computer. 
In fact the limitation to factors strictly not larger than one is only imposed 
for mathematical convenience. What is actually needed for convergence is 
that the functional be bounded by some value with probability one. A more 
relaxed condition on the constants may therefore be obtained by imposing 
|Pn£'max(^)^mail < ^ fo^ proccss I and jp^^^^^^ (//) F^axl < M for pro- 
cess II, pn being the probability of a tree with branchings, g"^^-^ (/) and 
fi'max i.^^) the maximum values of the couplings and -Fmax the maximum value 
of the initial condition. 

To test the method we have studied the time evolution of small and large 
Fourier modes in a plasma with two particle species of opposite charges, 
one light and the other heavy, with two types of initial Fourier distribution 
functions, namely 

i^/'H6,6,0) = C«e--l«^l^e-^»l«^l^ (32) 
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and 



i^f (6,6,0) 
i^!'^ (6,6,0) 



r7(2) -7i6i 



c(V^i«^i'e(fc-|6l') 



(33) 



= 40/3_ and the Cq s are chosen to fulfill condition (A). For the initial 



condition we have chosen = (0.01,0.01,0.01) and varied |^2|^ in the range 
3.10"^ to 7. Then, the time evolution is computed using the one-sided rep- 
resentation (process II). Some results are shown in Figs. 3 and 4, with time 
in units of j. Although it is known that on and without an external 
force there are no nontrivial steady states when both charges of opposite 
sign can move[TT] [12], one can see the relative stability of the Fourier modes 
for short times for the Gaussian initial condition F^^\ whereas for the F*^^^ 
initial condition one sees the appearance of growing Fourier modes that were 
not present in the initial density. Although the points were computed for the 
same set of final t's, in the plots we show the actual time t, obtained from 
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Figure 3: Time evolution of some Fourier modes for the F^^^ initial condition, 
T = Xt 
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Figure 4: Time evolution of some Fourier modes for the F^^^ initial condition, 
T^Xt 

Notice that to obtain a reasonable stability of the averaged functional one 
needs to compute many sample trees. The points shown in the figures were 
obtained with averages over 10^ or 2 x 10^ trees, depending on the evolution 
time. The reason for the need for a large number of sample trees arises from 
the fact that, for large times, most trees contribute a very small value to the 
average, the actual average arising from the contribution of a small number 
of them. This calls for the need to control the results by a large deviation 
analysis (see below). 

4 Conclusions 

1 - Mostly when localized solutions in Fourier space (or in configuration space 

for other stochastic representations) are desired, the method seems appro- 
priate. When a global calculation of the solution is desired, the stochastic 
representation method is probably not competitive with other current sim- 
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ulation methods. The computational example presented in Sect. 3, merely 
illustrative of the method, was obtained with modest computational means. 
Our purpose was mostly to test the stability of the results. To obtain good 
statistics and also to study the fluctuation spectrum of the process, many 
sample trees have to be used for each initial condition. However, because each 
tree is independent from the others and also because after the branching each 
mode evolves independently of the others, this algorithm is well suited for 
parallelization and distributed computing. In this sense the stochastically- 
based algorithms might also become competitive even for global calculations 
using parallel computing. In fact stochastic representations have already 
been found to be efficient for domain decomposition in parallel computing 

2 - The fluctuations around the mean in a branching process are typically 
very much non-Gaussian. Therefore a simple calculation of the standard 
deviation or other lower order momenta are not sufficient to check the relia- 
bility of the results. A large deviation analysis is recommended for numerical 
calculations using branching processes. Some general results on large devia- 
tions in branching processes are known [15j ^B] ^Tj [TS]. Of more practical 
importance are probably methods to estimate large deviation effects directly 
from the data. This may be done, for example, by the empirical construction 
of the deviation function. This is done by the empirical construction of the 
free energy and from it, by Legendre transform, the deviation function. For 
details we refer to [T^. Given a deviation function I (x), the probability of 
obtaining a value x for the empirical average of a sample of size n is 

where x means logarithmic equivalence. We have used the method described 
in [19] to check the reliability of the results. In the Fig. 5 we present the 
empirically obtained deviation function for a sample of 5 x 10^ trees. At first 
sight the regular behavior of / {y) around the mean, seen in the upper plot 
of Fig. 5, would seem to indicate that the distribution is Gaussian. However 
expanding a little more (in the lower plot) the domain of the variable x one 
sees the very non-Gaussian nature of the data. It means that, had we used 
a smaller sample, any empirical mean in the range 0.04 — 0.075 would have 
been likely. A rough lower bound on the size of the needed sample may be 
obtained from the inverse of the deviation function at the point where the 
behavior of / (x) changes. 
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Figure 5: The behavior of the deviation function for a sample of size 5 x 10^ 

3 - Stochastic representations of the solutions of deterministic equations 
may have some relevance for the study of the fluctuation spectrum. In the 
past, the fluctuation spectrum of charged fluids was studied either by the 
BBGKY hierarchy derived from the Liouville or Klimontovich equations, 
with some sort of closure approximation, or by direct approximations to the 
N-body partition function or by models of dressed test particles, etc. (see 
reviews in [20] [21] )• Alternatively, by linearizing the Vlasov equation about 
a stable solution and diagonalizing the Hamiltonian, a canonical partition 
function may be used to compute correlation functions [22] . 

As a model for charged fluids, the Vlasov equation is just a mean-field 
collisionless theory. Therefore, it is unlikely that, by itself, it will contain 
full information on the fluctuation spectrum. Kinetic and fluid equations 
are obtained from the full particle dynamics in the 6N-dimensional phase- 
space by a chain of reductions. Along the way, information on the actual 
nature of fluctuations and turbulence may have been lost. An accurate model 
of turbulence may exist at some intermediate (mesoscopic) level, but not 
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necessarily in the final mean-field equation. 

When a stochastic representation is constructed, one obtains a process 
for which the mean value is the solution of the mean-field equation. The 
process itself contains more information. This does not mean, of course, 
that the process is an accurate mesoscopic model of Nature, because we 
might be climbing up a path different from the one that led us down from 
the particle dynamics. Nevertheless, insofar as the stochastic representation 
is qualitatively unique and related to some reasonable iterative process, it 
provides a surrogate mesoscopic model from which fluctuations are easily 
computed. This is what we have referred elsewhere as the stochastic principle 
|lUj . At the minimum, one might say that the stochastic principle provides 
yet another closure procedure for the kinetic equations. 
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